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lilie low-le\'cl inteipretation of images provides constraints on 3D surface shape at multiple reso- 
lutions, but typically only at scattered locations over die visual field. Subsequent visual processing 
can be facilitated substantially if die scaticrbd shape constraints are immediately transformed into 
Visiblc-surfiice representations tliat unambigiiously specify surface shape at every image point. 'Hie 
required transformation is shown to lead tb an ill-posed surface reconstnicdon problem. A wcll- 
t^oscd variational principle formulation is obtained by invoking "controlled continuity,'* a physically 
flonicstrictive (generic) assumption about surfaces which is nonetheless strong enough to guarantee 
i^nique solutions. The variational principlcj which admits an appealing physical interpretation, is 
%ally discrelizcd by applying tlie finite el<|ment method to a piecewisc, finite element reprcsen- 
mion of surfaces. This forms the malheitiatical basis of a unified and general fra«ncwork for 
computing visible-surface representations, the computational framework unifies formal solutions 
to the key problems of (i) integrating multi^ale constraints on surface deptli and orientation from 
fliultiplc visual sources, (ii) interpolating difese scattered constraints into dense, piecewisc smooth 
siufaces, (iii) discovcritig surface depth and orientation discontinuities and allowing diem to restrict 
interpolation appropriately, and (iv) overcoi^ing the immense computational burden of fine resolu- 
tion surface reconstruction. An efficient surface reconstiuction algorithm is developed. It exploits 
multiresolution hierarchies of cooperative relaxation prcKcsscs and is suiUible for implementadon on 
massively parallel networks of simple, locally interconnected processors, 'lite algorithm is evaluated 
empirically in a diversity of applications. 
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I. Introduction 



Computing yisible^Surface Representations 



Over tliirty years ago, JJ. Gibson [1950] made the seminal conjecture tliat natural human perception 
amounts to the perception of visible surfhces. The explicit representation of visible surfaces, an 
intermediate goal of computational vision, has since attracted considerable interest. 

The computational framework offered in this paper addresses, in a unified way, certain visual 
information processing tasks involved in ttie representation of visible surfaces. Particular emphasis 
is placed on utilizing highly parallel, coopprative pnxressing to integrate surface shape information 
over multiple visual sources, to fuse it across a multiplicity of spatial resolutions, and to maintain the 
global consistency of the resulting distributed shape representations. The issues are first investigated 
fn terms of a surface reconstrucL mode] rooted in mathemaUcal physics, ll^is formal analysis is 
augmented by an empirical study of the resulting algorithms, which feature multiresolution iterative 
processing within hierarchical surface shabe representations. The approach is guided by current 
knowledge of how humans perceive visiblp surfaces, while applications in machine vision provide 
a testbed for die algorithms. 

The remainder of this introductory section examines the role of surface representations in early 
visual processing, outlines tlie key compi^tational problems that will be of primary concern, and 
reviews some relevant prior work. 



1.1. Early Visual Processing and Visible-Surface Representations 

Early vision comprises a set of processes ^hich specialize in recovering the physical properties of 
visible surfaces in a 3D scene from 2D imkes of the scene. They apply generic assumptions about 
the physical world and the imaging process to infer 3D surface shape constraints by interpreting 
specific image cues, such as stereoscopic disparity, motion, texture, contoure, and shading. These 
conceptually independent shape estimatioij processes fall into two broad categories. 



TTic first categorycomprises what are 
operate over multiple image frames of a 
arc stereopsis and structure from motion 
and [Ullman, 1983]). Stereopsis is driven 
simultaneously, but from different spatial 
involves frames taken from the same 
established across the frames, between 
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position, but at different times. If correspondences can be 

image features which originate from the same point on 

then the depth (i.e., 3D distance) to such points can 
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The second category of shape estimation processes involve computations on a single static 
frame. Perspective projection of 3D scene$ onto images imparts a systematic distortion to imaged 
surface properties such as shading, texture, and contours. A major part of tliis distortion can be 
attributed to the relative orientations of visible surfaces with respect to the viewer. In principle, 
it is possible to estimate surface orientation by measuring and interpreting such distortions in die 
image. This is the basis of practical approaches to recovering surface shape from shading, texture, 
and contours [Ikeuchi and Horn, 1981; Hoin and Brooks, 1985; Kender, 1980; Witkin, 1981; Brady 
and Yuille, 1984]. 

The combined output of the shape estimation pr(x:esses is best collected into intemiediate 
representations of die 3D shapes and coh figurations of visible surfaces, which we will refer to 
as visible-surface representations. Notable among proposed visible-surface representations are die 



Terzopoulos Compuiing Visible-Surface Representations 2 

depth and needle maps of Horn [1982], the intrinsic images of Barrow and Tcncnbaum [1978], and 
tlie 2|-1) sketch of Marr and Nishihara [1978]. For humans, the perception of visible surfaces is 
generally immediate, involuntary, and seems to precede (object) recognition. 'ITiis strongly suggests 
the existence of a visual process that autonomously computes visible-surface representations. Aside 
from the perceptual evidence, the availability of explicit visible-surface representations can also 
substantially facilitate subsequent surface analysis tasks in machine vision. 

Since early visual prcxessing provides relative surface shape estimates with respect to the 
viewer, it is most natural to define tlie basic shape primitives of visible-surface representations in a 
viewer-centered coordinate system. Moreover, the primitives should be computationally compatible 
with tlie local depth and orientation measurements (as well as discontinuities) that are provided by 
tlie various shape estimation processes. These criteria are satisfied by a particularly appealing class 
of local, piecewise shape primitives known as finite elements. 

A crucial realization is that shape estimates can be provided at multiple resolutions. Indeed, 
multircsolution spatial frequency channels have been identified psychophysically in tlie human 
visual system (e.g., [Braddick et at, 1978]). ITieir existence has influenced the design of early 
visual algorithms (e.g., [Marr, 1982]). In addition, machine vision research has demonstrated that 
multircsolution processing effectively bridges fine and coarse image structure, while it simultaneously 
increases computational efficiency (e.g., [Rosenfcld, 1984]). Hence, a multircsolution organization 
of visible-surface representations is most desirable [ferzopoulos, 1982, 1983a]. 

1.2. Key Problems of Visible-Surface Reconstruction 

The main topic of concern in this paper is the development of a visible-surface reconstruction process 
responsible for generating and dynamically maintaining visible surface representations. Whether 
the intention is to model human vision or to design competent artificial vision systems, this process 
must solve four key problems frcrzopoulos, 1983b, 1984]: (i) the constraint integration problem, 
(ii) tlic interpolation problem, (iii) the discontinuity problem, and (iv) die computational efficiency 
problem. We elaborate on each of these problems next. 

(i) The Constraint Integration Problem: Each specialized visual process may be thought of as a 
quasi-independent source of information partially constraining the shapes of visible surfaces. 
Tlic human visual system is reliable and robust because it integrates tlie various processes, 
enabling them to complement one another. I'he integration of multiple sources of information 
introduces redundancy, which is necessary not only to resolve potential ambiguities, but also 
to overcome the detrimental effects of noise and inaccuracies in the initial shape esdmates. 
'I'he constraint integration problem is fundamentally one of devising an effective means of 
integrating all available surface depth and orientation constraints (and discontinuities) within a 
cooperative visible-surface recomtruction process. 

(ii) The Interpolation Problem: It is widely accepted that initial descriptions of images ought to 
make explicit die occurrence and local 2D structure of image features diat are correlated to 
salient events on physical surfaces (markings, boundaries, etc.). This is the essence, for instance, 
of Marr's "primal sketch" representation of significant image irradiance changes (edges) [Marr, 
1982]. Generally, such salient features do not occur everywhere over tlie visual field. The 
initial representation of images as a sparse set of features implies tliat surface shape constraints 
generated by tlic specialized processes will also be scattered over a subset of image points. It 
is fascinating, however, that the human visual system systematically interprets visual stimuli 
such as sparse randoiri dot stereograms as coherent 3D surfaces [Julesz, 1971]. Indeed, tliese 
stereograms continue to elicit perceptions of dense surfaces, even when die density of dots 
carrying disparity information has been reduced until depth is unspecified over 98 percent of 
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the visible surfiicc area (see Fig. 1). It tiicrcforc appears tliat tlic surface reconstruction process 
is smoothly "filling in the gaps." This phenomenon has been the subject of some psychophysical 
investigation (e.g., [Collett, 1984]). The interpolation problem of visible-surface reconstruction 
challenges us to devise a scheme, consistent with human perception, for propagating shape 
infonnation into indeterminate regions (devoid of shape estimates) from places where it is 
available. 





Figure 1. A sparse random dot stereogram. Binocular fusion of this stereogram reveals a planar surface as a 
central, opaque, textured square suspended nearer in depth over a sunilarly textured background. Vivid depUi 
discontinuities separate the dense surfaces. 

(iii) The Discontinuity Problem: Visual discontinuities result from significant, spatially-localized 
changes in the physical world, particularly abrupt changes in surface structure. Botli depth 
and orientation discontinuities are perceptually relevant and provide vital boundary conditions 
for surface rcconstrucdon. Discontinuities in dcpdi occur at occluding contours, along which 
a surface in the scene occludes itself or another surface. Orientation discontinuities occur at 
creases or cusps of an otherwise continuous surface. In addition to the perception of coherent 
surfaces, random dot stereograms elicit vivid perceptions of surface discontinuities at abrupt 
disparity changes (see Fig. 1). The discontinuity problem amounts to (1) finding both dcpUi 
and orientation discontinuities in surfaces, and (2) dealing with their presence during visible- 
surface reconstruction; i.e., allowing discontinuities to limit tlic otlierwise smooth interpolation 
of shape constraints. 

(iv) The Computational Efficiency Problem: Visible-surface reconstaiction at die resolution of the 
image imposes an immense computational burden on both biological and artificial vision sys- 
tems. Nevertlieless, visible-surface representations must be computed quickly if diey arc to 
be of any practical value. It is generally accepted tliat to achieve die necessary performance, 
visual algoriUims and mechanisms must emphasize parallelism (liallard et al., 1983]; however, 
visible-surface reconstruction is compute bound to tlie point where the fundamental limitations 
of massively parallel mechanisms, particularly with respect to global interprocessor commu- 
nications, lead to severe inefficiencies. The computational efficiency problem is to develop a 
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visible-surface reconstruction process that not only exploits parallclisin, but also overcomes co- 
operative communication bottlenecks to compute visible-surface representations quickly, given 
suitable architectures. For the reasons outlined in the foregoing section and following our 
previous work [Icrzopoulos, 1982, 1983a], our solution to tliis problem hinges on the idea of 
multi resolution structuring of visual representations and associated cooperative processes. 

1.3. Prior Work 

There has been some prior work relating to tlie computation of visible-surface representations. 
Barrow and Tenenbaum [1979] describe an approach to rcconstiiicting smootli surfaces from noisy 
visual datii. This approach did not apply to general classes of surfaces, however, and the proposed 
relaxation algoritlims were not supported by a firm mathematical analysis. Nevertheless, Barrow 
and Fenenbaum's [1978] basic model of intrinsic images and much of the philosophy underlying 
tlicir computation seems appropriate, and it has influenced our approach. 

The interpolation problem is related to classical spline approximation. A number of well-known 
surfece approximation methods for scattered data are reviewed by Schumaker [1976]. Crimson 
[1983] employed one of Uiesc metliods for tlic continuous interpolation of visual surfaces from depth 
constraints; a minimization scheme involving a particular functional containing second derivatives 
(he referred to it as the "quadratic variation"). Brady and Horn [1983] observe tliat this functional 
is related to the bending energy of a tliin plate (a connection noted by Duchon [1977]), and tlie 
tliin plate model was developed further by Icrzopoulos [1983a] (see also [Blake, 1984]). 

Interestingly, thin plate interpolants have appeared in other areas, including the interpolation 
of aircraft wing deflections [Harder and Desmarais, 1972], interpolation of meteorological fields 
[VVahba and VVcndelberger, 1980], and the interpolation of digital terrain maps [Briggs, 1974; 
Bolondi et ah, 1976]. In this latter paper dicre is some concern for tlie presence of discontinuities 
(faults). 

Following Ullman [1979] and others. Crimson [1983] pursued "biologically feasible," parallel 
and iterative algoritlims for surface interpolation. A serious drawback of algorithms which satisfy 
these criteria is that they often converge cxcniciatingly slowly for problems of reasonable size, 
llie idea of multiresolution surfiice reconstruction exploiting multigrid relaxation mediods was 
shown to overcome diis problem while adhering to biological fcasibility [Terzopoulos, 1982, 1983a]. 
The multiresolution methodology yields efiicient algoriUims not only for the surface reconstrucdon 
problem but for other visual problems as well fFcrzopoulos, 1984]. 

In retrospect, akhough progress has been made, a satisfactory compuUitional theory of visible- 
surface representations has been elusive. This is largely a consequence of the significant technical 
obstacles encountered in devising fomial solutions to all four key problems of visible- surface 
reconstruction within a unified computational framework. The difliculty of tlie task appears to have 
evoked some skepticism as to the actual computability (hence, even tlie usefulness) of intrinsic 
surface representations [Witkin and' Tencnbaimi, 1983]. Based on the theoretical generality of 
our approach and the accompanying empirical' results, however, we believe such skepticism to be 
premature. 



2. Mathematical Analysis of Visible-Surface Reconstruction 

Let the true distance from the viewer to visible surfaces be given by the function Z(x,y), where 
X and y arc the image coordinates. Fow-level visual processes generate a set of noise corrupted 
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surface shape estimates (i.e., constraints) {q} whicli can be expressed in the abstract notation 

c^ = Jli\Z(x,y)]-^U, (1) 

where Ci denote measurcnncnt functionals of Z(:r,'i/) and ({ denote associated measureinent errors. 
Stated simply, the visible-surface reconstruction problem is to reconstruct, as faithfully as possible, 
tlic depth function Z(x,y) from the available constraints {c,}. 

2.1. 1 he Ill-Posed Nature of the PVobleni and Rcgularization 

The problem is made nontrivial by the nature of the constraints. First, constraints are contributed not 
by one, but by multiple specialized early visual processes. Hence, slightly inconsistent measurements 
provided by different processes that happen to coincide will locally overdetcrmine surface shape. 
Second, constraints are not dense, but scattered sparsely over tlic visual field. Therefore, while 
diey may restrict surfiice shape locally, tlicy do not determine it uniquely everywhere; Uierc remain 
very many feasible surfaces that are consistent with the constraints. Third, the measurements are 
subject to errors and noise. High spatial frequency additive noise, regardless how small its (RMS) 
amplitude, can locally perturb the surfoce (orientation) radically. 

In view of the above three considerations, we cannot conclude in general that the solution 
will exist, nor that it will be unique, nor that it will be stable with respect to measurement 
errors. Mathematical problems for which the existence, uniqueness, or stability of solutions cannot 
be guaranteed a priori are said to be ill-posed [Tikhonov and Arscnin, 1977]. Visible-surface 
reconstruction can tlius be characterized as a fundamentally ill-posed problem. 

Ill-posed reconstruction (or inverse) problems are the rule radier than die exception in early 
vision [Poggio and Torre, 1984J. Ill-posed problems cannot be solved in general, without imposing 
some additional restrictions on possible soludons. This is die basis of a number of systematic 
approaches, notably the rcgularization methods introduced by Tikhonov and others (see [Tikhonov 
and Arscnin, 1977] and references dierein). Duda and Hart [1973, Sec. 7.4] mendon a basic forni 
of rcgularization (essentially spatial smoothing) for combating die effects of noise in images. A 
more sophisticated class of rcgularizadon methods is discussed in the context of low-level vision by 
Poggio and Ibrre [1984]. 

Through rcgularization, ill-posed problems can be solved by reformulating them as mriational 
principles that are effectively computable. Unlike the original problems, die variational principle 
formulations are well-posed; i.e., it is possible to guarantee the existence, uniqueness, and stability 
of dicir solutions under nonrestrictive conditions. Reformulation proceeds with die introduction of 
suitable slabilizing funclionals, notably the class of stabilizers proposed by Tikhonov and Arsenin 
[1977, pp. 69-70]. These stabilizers can be interpreted as spline finicdonals that impose smoothness 
assumptions on the admissible solutions (by restricdng them to Sobolev spaces of smooth functions).-^ 
Pragmatically then, this type of rcgularization is essentially equivalent to optimal approximation by 
generalized splines [lerzopoidos, 1985a]. We pursue die generalized spline approximation point of 
view, since splines are familiar and since they suggest helpful physical interpretations. 

2.2. A Variational Principle 

The abstract dieory of optimal spline approximation is well-developed and a close connection has 
been established with variational principles involving the constrained minimization of (semi-) norms 
in (semi-) Hilbcrt function spaces [Laurent, 1972]. Let ^ be a linear space of smooth functions and 
let S{v) be a functional defined on ^ which measures the (lack of) smoodiness of a function in 



^ Generic smoothness assumptions are generally the weakest (least committal) assumptions that one can make 
about feasible solutions and still obtain well-posed formulations. 
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)i. Furthermore, let /' be a functional on 'H which provides a measure of tlie discrepancy between 
the function and the given constraints. Consider the following variational principle: 

VP: Find u e M such that 

^W = inf^W, (2) 

where the energy junctional 

£{v) = S(v) + P{v). (3) 

This variational principle will serve as a formal statement of Ihe visible-surface reconstruction 
problem: The best reconstruction of the deptli flmction Z(x,y) from die available constraints will 
be given by the solution u(x,y), tlie smoothest function in die admissible space ^ which is most 
compatible wiUi die constraints. 

Before proceeding to specify die smoothness functional S{v) and die penally Junctional P(v), 
it should be noted diat, if the solution exists, it satisfies die necessary condition for die minimum 
given by the vanishing of die first variation 6, 

6£{u) = 6S(u)+6P(u)=0, (4), 

which expresses the so called Euler- Lagrange equations. 

2.3. Generalized Spline Functionals 

For an appropriate smoodiness functional 5(w), we turn to die muUidimensional sphnes studied 
by Duchon [1977] and Meinguet [1979], generalizations of die classical univariate splines [Ahlberg, 
et al, 1967]. The subclass of {2D) surface splines relevant to our problem can be characterized as 
members of a suitable space of admissible functions v(x,y) which minimize die functional 

|<=/lS(7)(^-)'-- 

ITie posidve integer m dictates die order of the pardal derivatives that occur in die funcUonal, 
which in turn determines die order of continuity possessed by die admissible functions. The Euler- 
Lagrange equadon satisfied by the minimizing function u(x,y) is an iterated version of Laplace's 
equation: (-l)''"A^"u = 0, where Aw = u^^ + Uyy is die Laplacian of u. 

Low order surface splines have interesting physical interpretations involving equilibria of elastic 
bodies. Two special cases arc of interest. For m = 1 die functional reduces to 



Ml = j J K' + "') «'^^^2/, (6) 



which is proportional to the small deflection energy of a membrane (e.g., rubber sheet), while for 
m = 2, 

•2 +2^2^ + uJJ dxdy, (7) 



//fe 



is proportional to the small deflection bending energy of a diin plate (with zero Poisson ratio) 
[Courant and Hiibert, 1953]. Duchon [1977] refers to the minimizers of \v\l as tiiin plate splines. 
Since diin plate splines are the natural 2D analogs of cubic splines, \v\l finds frequent usage in 
surface inteipolation problems [Schumaker, 1976]. In particular, it has been employed for visual 
surface interpolation [Crimson, 1983; Terzopoulos, 1983a]. 

The physical interpretations make it clear that membrane splines offer a lower order of 
continuity dian thin plate splines. Since die physical forces in die membrane are due primarily 
to its surftice tension, it generates minimal area surfaces. Although minimal area surfaces are 
continuous, they need not have continuous first partial derivatives; i.e., they are C^ surfaces. For 
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insUincc, a sharp corner would result readily if an idealized physical membrane were subjected to 
tlie deflecting force of a knife edge. In contrast, the restoring forces in a physical thin plate are 
due primarily to its flcxural rigidity. A thin plate would not crease when deflected by a knife edge. 
Thin plate splines tJierefore maintain continuity as well as continuous first partial derivatives; i.e., 
they generate C^ surfaces. 

2.4. Controlled Continuity and the Thin Plate Surface Under Tension 

Generic smoothness assumptions are justified in pursuing a regularization approach to die visible- 
surface reconstruction problem, inasmuch as the coherence of matter tends to give rise to smoothly 
varying surfaces relative to the viewing distance, over some range of scales; however, smoothness 
assumptions clearly do not hold arbitrarily across surface discontinuities, some of which persist 
across all scales. This introduces significant complications for classical spline approximation or 
regularization methods; the continuity of spline fimctionals (or stabilizers) must be controlled at 
discontinuities in order to preserve them. 

A stabilizer providing the necessary local continuity control can be realized as a weighted 
combination of generalized spline fimctionals of more than one order m (Terzopoulos, 1985a]. We 
propose the following smoothness functional: 

Spr(v] = lj J p(^>2/){^(^,y)(^L + 2t;2^ + ^yj/) + [i - r(x,y)](vl + vl)}dxdy, (8) 

where denotes the image domain, and p(x,y) and T{x,y) are real-valued weighting functions 
whose range is [0, 1]. This controlled-continuity stabilizer is a weighted convex combination of the 
thin plate spline functional \v\l and membrane spline functional \v\l integrands. The associated 
Euler-Lagrange equation is 

dx^ *^''=^='' ■^ d^ '^^''"^' '^' a^2 i^^^yy) - ^ (^^x) - ^ (»?%) = 0' (9) 

where //(x, y) = p{x, y)T{x, y) and r/(a;, y) = p{x, y)[l-T(x, «/)], with natural (i.e., free) boundary 
conditions. The flinctional Spriv] can be thought of as a thin plate surface under tension, where 
p{x,y) is a spatially varying "rigidity" and [1 - T(x,y)] is the spatially varying "surface tension." 
It generalizes the unidimensional splines under tension of Schwcikcrt (see [Ahlberg, et ah, 1967]). 

The local continuity properties of the thin plate surface under tension functional can be 
controlled at any point (x,y) G n by specifying tlie values of the continuity control functions 
p{x,y) and T{x,y) at tliat point. As r approaches 1 the fiinctional tends to a thin plate spline (a 
C^ surface) whereas towards tlic other extreme, 0, the functional tends to a membrane spline (a 
C^ surface) witli intermediate values characterizing a hybrid C^ surface tliat blends tlie properties 
of both constituent splines, p determines the overall potency of die smoothness functional. 

Reconstnicted surfaces must be able to faithfully preserve known depth and orientation 
discontinuities, while not introducing spurious discontinuities at other locations. This can be 
accomplished if (i) away from known dcptliand orientation discontinuities, die reconstructed 
surface possesses (at least) the C^ smoothness of a tliin plate, maintaining both continuity and 
continuous derivatives, (ii) at known orientation discontinuities, it exhibits just die C^ smoothness of 
a membrane, maintaining continuity only, and (iii) at known depth discontinuities, the smoothness 
functional is deactivated so that the reconstructed surface is free to "fracture" locally. Hence, 
Spt(v) will be manipulated as follows: At all non-discontinuity points (x,?/), p(x,y) and T(x,y) 
should be nonzero. At orientation discontinuity points, T(x,y) is set to zero. At depth discontinuity 
points, p{x,y) is set to zero. Mechanisms for aiitomatically detecting discontinuities by computing 
continuity control functions optimally according to local criteria are considered in a subsequent 
section. 
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2.5. Penalty Functionals 

Assuming independently distributed measurement errors c,: witii zero means and variances a?, tlic 
optimal measure of incompatibility is a weighted Huclidean norm of tlie discrepancy between the 
admissible function and the data q: 

^w=ii;".(AH-'^.f, (10) 

where the ai are nonnegative real-valued weights (ideally «,■ is inversely proportional to (t?; - 
i.e., ai = 1/Aa?) [Kimeldorf and Wahba, 1970]. This penalty functional can also be employed 
(suboptimally) when the above assumptions do not hold strictly. 

Appropriate measurement functionals Li for surface reconstruction may be synthesized from 
generalized A;^-order derivatives: 

^•M=5^L' > = M,....fc. (U) 

k = yields simple evaluation functionals Cilv(x,y)] — v(xi,yi), which will be employed to 
model the local depth constraints 

Ci = v(xi,yi) + Ci = 'i(x,,y,)- (12) 

The components of tlie local surface normal n(xi,yi) — \v^[xi,yi)^Vy[xi^yi\-\\, which deter- 
mine local surface orientation, can be handled by the first order (/c = 1) derivative functionals 
Li\x){x,y)\ - Vx{xi,yi) and Cilv(x,y)] = Vy(xi,yi) and yield analogous expressions for the local 
orientation constraints: 

Ci =Vy{Xi,yi) +€i =q{xi,y,)- 

Other potentially relevant functionals such as directional derivatives can be accommodated straight- 
forwardly with the above notation. 

It is convenient to separate tlic various constraints into three sets; the set i e D of image 
points at which depth constraints cf(x.,i/,) occur, and the sets i e P and « G Q at which orientation 
constraints p(j;.^y.) and ^(i.,^.) occur respectively. The penalty functional can then be expressed as 
a sum of three components 

ieD ieP ieQ 

(14) 
where the a,- parameters are now distinguished as aji, ctp., and aq.. 

2.6. A Physical Model for Visible-Surface Reconstruction 

llic variational principle formulation of the surfocc reconstmction problem has an appealing 
physical interpretation which is illustrated in Fig. 2. The thin plate surface under tension may 
be visualized as an elastic surface, planar in its natural state, whose clastic bending energy Spt(v) 
stabilizes surface shape so tliat it varies smoothly in between constraints (but not at discontinuities). 
Constraints deflect the surfece according the penalty functional P(v), which can be intcipreted as 
tlic total stretching energy of a set of ideal springs attached to tlic constraints. The left part of 
tlie figure shows the elastic surfoce whose deflection u{x,y) at equilibrium is determined by an 
infrastructure of scattered dcptli constraints. The local depth estimate is encoded as the vertical 
height of the constraint and the tightness of each constraint is controlled by associated spring 
stiffness ocdi. The right part of tlic figure illustrates an orientation constraint coercing tlie local 
surface normal. The spring stiffness is determined by the constraint parameters (Xp- and rv^-. 
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orientation constraint 



surface normal 




Figure 2. The physical model. Thin plate surface under tension and depth constraints (left). Local influence 
of an orientation constraint (right). 



2.7. Existence, Uniqueness, and Stability of the Solution 



Existence, uniqueness, and stability of the solution u(x, y) to the variational principle VP are 
guaranteed when Spriv) = Spriv) + P(v) is a norm in tlic admissible space ^. Unfortunately, 
generalized spline functionals \v\1^ are a priori only semi-norms (of a particular class of Sobolev 
spaces). The null spaces .V of functions tliat map to zero under tlie semi-norms are simply the 
(M = C"^^) dimensional) spaces of all polynomials over 9?^ of degree less than or equal to m- 1 
[iMeinguet, 1979]. The penalty functional P(v) can force Spriv) to be a norm, however, if it at 
least constrains .A/ to a unique polynomial. A possible set of conditions for this to occur is tliat the 
Ci include evaluation functionals at an >/-unisolvent set of points (i.e., a set of M points which 
define a unique polynomial in the null space of tlie smoothness flmctional). In particular, since the 
maximum order of generalized splines in the stabihzer Spriv) is m - 2, its null space is the space 
of linear polynomials. Thus, tlie following proposition can be proven [Terzopoulos, 1984]: 

Proposition. The solution u(x, y) will exist, be unique, and stable given any one of the following 
minimal conditions 

(i) three noncolinear depth constraints, 

(ii) two depth constraints as well as a single p or q constraint, 

(Hi) a single depth constraint as well as a single.p and a single q constraint, 

(iv) a single p and a single q constraint with the "center of gravity" of the surface fixed. 

These minimal conditions will hold in practice, due to the large number of constraints typically 
available from early shape estimation processes (tlie fixed center of gravity condition can be imposed 
when necessary). Consequently, die visible-surface reconstruction problem may be considered well- 
posed, hence effectively computable in general. 

Satisfying the conditions for a well-posed problem essendally guarantees that a unique state of 
stable equilibrium will exist lor the plate/spring system (the minimal energy sUite Spriu)). In this 
context, the controlled continuity assumption about surfaces, as embodied by the thin plate surface 
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under tension model, is physically nonrcstrictive but nonetheless powerful enough to guarantee tlie 
existence of unique solutions to tJie variational principle. 



3. Discretization 



It is extremely difficult, if not impossible, to obtain an analytic solution to the variational principle 
due to tlic irregular occurrence and geometry of constraints and discontinuities. For our purposes, 
the only viable approach is to convert tlie continuous surface reconstruction problem to an equivalent 
discrete problem whose solution can be computed numerically. To diis end, finite elements make 
ideal local surface shape primitives for use in visible-surface representations [Terzopoulos, 1982, 
1983a]. The finite element mctliod [Strang and Fix, 1973] is a general, powerful, and matliematically 
rigorous approximation technique which guides tlie selection of appropriate elements and governs 
their interactions according to the nature of the variational principle. 

The finite element method offers substantial flexibility in discretizing domains with irregular 
shaped boundaries. Although the use of irregularly shaped elements to discretize such domains 
may not present a feasibility problem with regard to distributed biological mechanisms, it makes 
nontrivial the mapping of elemental computations onto regularly interconnected processing networks 
typically provided by VLSI technology. In this paper we restrict ourselves to regular finite elements 
in order to facilitate such mappings. Since the goal is to obtain a particularly fine discretization, at 
the resolution of the image, tlie restriction to fine regular elements will not jeopardize our abifity 
to accommodate tlie irregular occurrence of constraints or discontinuities. 



3.1. The Discrete Equations 



The domain (1 is tessellated into square element subdomains with sides of length h. Nodes are 
located at element corners and shared by adjacent elements. This results in a planar and uniform 
square grid of nodes tliat is ideally suited to VLSI implementation. Ihc nodes are naturally 
indexed by (i,j) for i = i,...,N^ and j = l,...,Ny, where N^ and Ny are tlie number of 
nodes along the x and y axis respectively of the (rectangular) domain Q. The total number of nodes 
is N - Nx X Ny. The reconstructed surface is represented by an assembly of (nonconforming) 
finite elements, each of which is a six-point (hill) quadratic intcrpolant defined locally widiin its 
particular subdomain. The unknown displacement (surface dcptli), at node (ij) is denoted by tlie 
variable v^^ = v(ih,jh). Taken together, the displacement variables are denoted by the vector 

v^ G 9ff^. Once tliis vector is determined by solving a discrete version of the variational principle, 
the local interpolants are known exactly and, consequently, they explicitly represent deplli and 
orientation everywhere over the surface. 



The proposed square, quadratic element leads to the following 0(/t^) fomuilas for the required 
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partial derivatives at an arbitrary node (z,j) ficrzopoulos, 1983a]: 

.h 






yy 



v„ 



xy 






(15) 






Note tliat the formulas arc finite difference expressions. Their appearance is due to the uniformity 
and low order of the element and dicir relative simplicity will facilitate the calculations substantially. 

Substituting the above expressions with the constant approximations p[x, y) => p^ ■ and 

t(x, y) => r^. into (8), and noting diat the area of each element is /i^, we obtain tlie discrete 
ftinctional 



s>'') = lEpt 



tj 






+ 



hj 






(16) 



Although by no means a necessity, it is both natural (in view of common image discretization) 
and convenient to assume that tlie constraints coincide with nodes (i,j) of the grid. Hence, to 
obtain a discrete expression for P(v), we collect the nodes at which the various constraints occur 
into three sets; the set (i,j) G D at which depth constraints d'',,- occur, and the sets [ij) G P 

and (i,j) G Q at which orientation constraints pfy and q^j occur. Using symmetric difference 
approximations for the partial derivatives in (13), the discrete penalty functional may be written in 
terms of the nodal variables as 

{i,j)eD 



2 ^ 









(17) 



(»J)£Q 

The energy-minimizing vector of nodal displacements u^ satisfies the equilibrium condition 

V^/Uu'^) = VSJUu'^) + V;''^(u'^) - 0, (18) 

where V is the gradient operator. Since the discrete functional Sp^iu^^) is a quadratic fomi in the 

u^j, the above equation defines a linear system of simultaneous equations that are satisfied by u'^. 
'I'he discrete problem amounts to solving these nodal equations. 
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To progress towards explicit expressions for the nodal equations, we first determine the partial 
derivatives of Spr{^^) and P^(w^) with respect to an arbitrary nodal variable uj*y. Letting 

(19) 



we obtain 



<3 = <3<3l^ ^"d r7,ti = P?,i(l-^^)' 



<3 



J-1 
1 



+ 



Ki - 2"?-i,i - 2u?,y_i + 2u?_i,,._i) /^tl, 
-2u?^.i,y + 2u?,,. + 2uti,y_i - 2u?,,._i) nl 
-2ut-+i + 'i^luHi + 2u?,y - 2uti,,.) M?--i„ 
2u?+i,y-,i - 2u^,.,i - 2u?+i,,. + 2u?,,.) ^ll^ 

<3 - 2"'i-i + <y-2) /^!;i-i 
-2u?,y^.i + 4u?,.-2u?,,._l)M!;y 



+ 



the discrete version of (9). Next, for (i,j) e Dn PnQ, 






'^^ij 



4.(^^A:l2i(,,h ..h \ ^Pi-l ,j h 

^[ 4/i2 V"»'^. "«+2,yj+ 2/z P»+i'^' 
^[ 4fi2 r»>^ "'.J~2J 2/1 "^^'^-i^ 



(20) 



(21) 



+ 



The above expressions specify tlie nodal equations implicitly. Each constituent tenn in (round) 
parentlicses can be represented graphically as a basic computational molecule. Computational 
molecules will be interpreted both as spatial representations of tlie nonzero coefficients in tJie nodal 
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equations, and as local computations involving multiplications and additions of specific proximal 
nodal variables. The former interpretation will facilitate the construction of individual nodal 
equations given some local structure of constraints and discontinuities, while the latter will lead 
directly to local iterative algorithms for solving the resulting simultaneous linear system. 



3.2.1. Basic Molecules 




1 w-2 v-ri 




-2V-(r4l-r-2 




1 M-2 w 1 










©-0 0-© 





Figure 3. Plate molecules (a) and membrane molecules (b). 
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Mq. (20) is a convex combination of two components; tiic first stemming from the thin plate 
energy functional, the second, from the membrane energy functional. Hach constituent term yields 
a basic computiUional molecule (sec Fig. 3). a set of linked atoms indicated by circles. The central 
node {i,j) is indicated by a double circle in each molecule. Fig. 3(a) illustrates the ten plate 
molecules obtained from die terms of the first component, while Fig. 3(b) shows the four membrane 
molecules obtained from tlie tenns of the second component. F^ch atom contains the coeflRcient of 
the associated nodal variable (aside from the fi^j and rj^j factors). 




3^ d^ 



4h^ 




-1 



a 







-1 



2h P-^'^ 



2/1 P»+i-' 







a * 




2h "^''-^ 



2h ^*''+^ 



Figure 4. Depth constraint molecule (a) and orienlaUon constraint molecules (b). 

Similarly, the depth constraint tenn in (21) can be represented by the depth constraint molecule 
shown in Fig. 4(a). Associated with it is the factor nJ^A^^, which is indicated underneath the 
molecule. The remaining orientation constraint tenns of (21) are represented by the orientation 
constraint molecules and associated factors shown in Fig. A{b). 
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The formation of nodal equations within continuous regions can be visualized as a process of 
molecular summation. During molecular summation, tlie basic molecules combine at the central 
node, coincident atoms summing together. 

When (i, j) is an interior node, away from constraints and discontinuities, p^- = r/\ = 1, and 
only the plate component of (20) contributes to the expression for the partial derivative. Hence, 
the equilibrium condition (18) reduces to the nodal equation 

^ (u?-i,i + u!Vi,,-fu^,._i + u?,^.+i) 
+ ^ ("ti,i-i + "iVi,y-i + "?-i,y-fi + »!Vi,y+i) (22) 



« 20 h 

= — u- 



+ p ("t2,y + U?+2,y + U?j.-2 + U?,y+2) 



This equation can be represented by the composite nodal molecule illustrated in Fig. 5(a), which 
results from die summation of the plate molecules in Fig. 3. 





Figure 5. Interior node molecules, (a) Away from discondnuities. (b) At interior orientation discontinuities. 

Note tliat the computational molecule for the center of the region is a factor of h"^ (due to die 
elemental area) times an order 0(/«^) finite difference approximation for tlie bihamionic operator 
[Abramowitz and Stegun, 1965, p. 885J, the Euler-Lagrange equation associated with the tliin plate 
spline. This is an expected consequence of tlie particular element employed which yielded finite 
difference approximations for the second partial derivatives of v^. 

If node (i,j) is a depth constraint, the first term in (21) takes part in the nodal equation. The 
effect can be represented as a summation of the deptli constraint molecule and associated constraint 
factor with the nodal molecule for {i,j) shown in Fig. 5(a). 

Similarly, if (i - l,j) or [i + IJ) arc p constraints, or (i,j - 1) or {i,j + 1) arc q constraints, 
the other terms in (21) participate in the nodal equation. Again tliis can be represented as the 
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summation of compuUUional molecules. Specifically, the upper left molecule in Fig. 4(b) sums with 
the nodal molecule if (i - ij) G P, tlie upper right only if (i + 1, j) G P, tlie lower left only if 
{i,j - i) eQ, and the lower right only if {i,j + 1) 6 Q. 

3.2.3. Molecular Inhibition at Discontinuities 

If (i,y) is a discontinuity, either ^i^- or 77-^^ or both may be zero, thus nullifying the summation 
of specific molecules. This crucial influence of the discrete continuity control flinctions near known 
discontinuities will be referred to as molecular inhibition. It was convenient for expressing (20) and 
(21) to discretize u(x,y), p(x,y), and r(a;,y) over die same set of nodes. Although orientation 
discontinuities can be situated at these nodes (since u^j is defined at an orientation discontinuity), 

it is better to position depth discontinuities on tlie links half way between nodes (since uj* • is 
undefined at a depth discontinuity). Discretizing T(x,y) on links does not present a problem in 
practice. As a general rule, a discontinuity may inhibit a molecule only if it coincides with a 
constituent atom or link. 

First consider orientation discontinuities. At an orientation discontinuity, r/* • = and only 
the second component of (20) contributes to the nodal equation. In effect, the plate molecules are 
inhibited and replaced by the membrane molecules of Fig. 3. At an interior orientation discontinuity 
(z, j), away from depth discontinuities, all four membrane molecules superpose to yield the nodal 
molecule shown in Fig. 5(b), which represents tlie nodal equation 

4u?,y - u?_i,,. - u?.,i,y - u?,y_i - u?,,. ,1 = 0. (23) 

llie equation will be recognized as -h? times a standard finite difference equation for the Laplacian 
[Abramowitz and Stegun, 1965, p. 885]. It too appears as a consequence of the Euler-Lagrange 
equation associated with the membrane spline. 

Since an orientation constrain cannot meaningfully coincide with an orientation discontinuity, 
orientation discontinuity nodes inhibit orientation constraint molecules. On the other hand, depth 
constraint molecules are not inhibited by orientation discontinuities since it is perfectly reasonable 
to locally constrain a membrane spline in depth. 

Because smoothness constraints are unsuitable at a depth discontinuity node (i,j) (i.e., /?t,y = 
0), a nodal equation cannot involve nodal variables separated by or coinciding with a depth 
discontinuity. Consequently, depth discontinuities inhibit all of tlie basic computational molecules. 
Fig. 6(a) illustrates examples (disregarding constraints) of nodal molecules for boundary nodes 
(marked as double circles) which are near depth discontinuity nodes (marked by X's). Examples of 
nodal molecules at boundary orientation discontinuities (double circles) next to depth discontinuities 
(X's) are shown in Fig. 6(b). 



4. Detection and Localization of Surface Discontinuities 



An important feature of our framework for computing visible-surface representations is the uni- 
fonn treatment of constraints and discontinuities, essentially as localized and independent surface 
shape primitives. This facilitates the parallel integration of discontinuity infonnation, along with 
shape constraints, over the various early shape estimation processes. It is convenient to think of 
discontinuity information as being collected into a discontinuity map which is in registration with 
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Figure 6. Molecular inhibition at discontinuities, (a) Nodal molecules at boundary nodes (double circles) near 
depth discontinuity nodes (X"s). (b) Nodal molecules at boundary orientation discontinuities (double circles) 
next to depth discontinuities (X's). 
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the reconstructed surface. Technically, the map comprises the nodal variables {p^A and {r^M 
representing the discrete continuity control functions. 

Any early visual process can participate in initializing tlie discontinuity map according to its 
own local hypotheses about the occurrence of discontinuities. In general, tliis prior discontinuity 
information will be partly incomplete and inconsistent, since it derives from narrowly specialized 2D 
image analysis. In evolving a globally consistent surface, the visible-surface reconstruction process 
performs a cnicial task: it brings the prior discontinuity information into consonance with the 3D 
shape constraints collected from all tJic early processes. This raises the problem of detecting and 
localizing both depth and orientation discontinuities during reconstruction. The current section 
investigates this problem for the (impoverished) case in which no prior discontinuity information 
is available. We first propose a straightforward scheme which exploits the regularizing properties 
of the surface model, then a more sophisticated approach that extends the variational principle to 
optimally estimate discontinuities according to generic expectations about their local structure. 

4.1. Regularization Based Discontinuity Detection 

From one perspective, surface discontinuity detection shares much in common witli traditional 
approaches to image intensity edge detection. In particular, it is possible to detect discontinuities 
by applying thresholded local differencing operations to the reconstructed surface which, like the 
image, is a regularly sampled function. Because they arc easily corrupted by image noise, however, 
local edge operators such as Laplacians perform poorly [Roscnfeld and Kak, 1982] without a 
smoothing prefilter, say a Gaussian fMarr and Hildreth, 1980]. Interestingly, the thin plate surface 
under tension performs the necessary smoothing on the sparse and noisy shape constraints (standard 
low-pass filters such as Gaussians are inapplicable to sparse data). This regularizing effect permits 
the reliable computation of numerical derivatives for detecting discontinuities [Bakhvalov, 1977, 
Sec. 5.4; Poggio and Torre, 1984; Terzopoulos, 1985a]. In addition to exploiting the regularizing 
effect of tlie diin plate surface under tension, tlie discontinuity detection scheme described next 
is easily accommodated within the distributed computational stRicture of our framework, and it 
permits relevant criteria such as psychophysically measured limits on stereoflision to impact on 
discontinuity detection. 

Consider the random dot stereogram in Fig. 7 which depicts a set of planar surfaces stacked in 
depdi. Fig. 8 shows a single continuous surface generated by the surface reconstruction algorithm 
from sparse stereoscopic disparities provided the Marr-Poggio-Grimson (MPG) stereo algorithm 
[Grimson, 1985]. Fig. 9 dramatizes a portion of the reconstructed surface in cross section as it passes 
across a depth discontinuity. The C^ surface overshoots constraints near the discontinuity because 
its smoothness conflicts with the sudden change in depth. The surface is clearly inappropriate as a 
final solution near deptli discontinuities, but the local incompatibility can signal tlie occurrence of 
tliese discontinuities. 

Opposing bending moments are imparted to the surface by the constraints on eitlier side of 
the discontinuity, llie surface inflection (see Fig. 9), where tlie bending moment undergoes a sign 
change, localizes die deptli discontinuity. For a thin plate spline u(x^y), the bending moment per 
unit length parallel to the x-z plane is proportional to -Uxt, while its counterpart parallel to the 
y-z plane is proportional to -Uyy [Szilard, 1974]. The sum of orthogonal bending moments gives 
die total moment M = -{uxx+Uyy) = -Aw, the negative Faplacian of the deflection function. It 
can be computed readily at a node (?', j) of the discrete surface using the standard approximation: 

M,i = -^ {ul,j + u?,i, • + uIj_ 1 + u?,,. , 1 - 4u?,y) . (24) 
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Figure 7. Synthesized random dot stereogram. When fused, the stereogram depicts four planar surfaces 
stacked one atop the other in depth. 



The zero crossings of M for the reconstructed surface in Fig. 8 are shown on the left in 
Fig. 10 as black contours. Most of these correspond to weak inflections due to slight ripples in tlie 
reconstructed surface. A measure of significance is dierefore needed to detect tine discontinuities 
while weeding out spurious, weak inflections. The magnitude of the local depth gradient (surface 
tilt) is a suitable significance measure for depth discontinuities. Hence, an inflection point will 

be considered significant if (7 = \Vu\ =■ Jul + uj exceeds a limit td (it is more efficient to 

use the square of this expression u^ + Uy, or even l^xl + |%|). Employing the usual discrete 
approximations, we obtain 



^'■' = iJ^ [("'+!.; - "U,y+ (<y+i - '''y-i)'] 



(25) 



The right half of Fig. 10 shows the significant inflection points where Gij > td with td = 1. 
Adding these significant points to the discontinuity map (by setting tlic associated pl^j to zero) 
fractures the continuous surface to yield as a solution tlie reconstructed stack of surfaces shown in 
Fig. 11. 

ITie limit td must be large enough so that weak inflection points are rejected as possible 
discontinuities, while not so large as to miss many true depth discontinuities. A possible criterion 
for choosing td in applications to stereopsis of opaque surfaces is suggested by Panum's limiting 
case: i.e., when a surface is tilted so much from tlie viewer tliat it begins to occlude itself from one 
eye, causing stereopsis to fail. Himian stereofusion limits have been measured psychophysically. 
Using pairs of points at different orientations, Burt and Julesz [1980J measured a roughly isotropic 
disparity gradient limit of approximately 1 between fusion and diplopia. Interestingly, this is only 
half the Panum limit. 

It is not inconsistent with these findings to use the disparity gradient limit td to detect 
significant depth discontinuities in conjunction with the isotropic bending moments M,;j to localize 
tliese discontinuities. The required local support computations can be performed in parallel at each 
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Figure 8. Single surface reconstructed from the stereogram of Fig. 7. 



grid point over the surface. Analogously, significant orientation discontinuities may be detected 
when the magnitude of tlie bending moment \Mij\ of the surface exceeds a limit to (points of 
high curvature), and they may be localized at relative extrcma of the bending moment (positions 
of locally highest curvature). The sign of a bending moment extremum indicates the sense of the 
orientation discontinuity; negative signals a concave crease, and positive, a convex crease. Curvature 
peaks were also employed in a scheme for detecting surface orientation discontinuities proposed by 
Langridge [1984]. 

4.2, Discontinuity Detection by Variational Continuity Control 

On the one hand, experimentation on natural data with the regularized approach to discontinuity 
detection demonstrates the feasibility of discovering many of the more significant discontinuities 
during surface reconstruction (results are presented later). On the other hand, certain inherent 
inadequacies of this simple scheme can often lead to poor surfece reconstaictions. The shortcoming 
are due to a basic conflict caused by smoothing. While rcgularization eliminates noise, making 
reasonable estimation of surface derivatives possible in continuous regions, it tends to obscure 
discontinuities [Terzopoulos, 1985a]. It can result in poor detectability and localization of the more 
subtle discontinuities, a common problem with smoothing edge operators in general [Leclerc and 
Zucker, 1984], 
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^ significant inflection 



V insignificant inflections 




Figure 9. Cross-section of a reconstructed surface across a depth discontinuity. The significant and insignificant 
inflections of the surface are indicated. 



The problem can be resolved by exploiting more fiilly the controlled continuity model to 
preserve surface discontinuities and, moreover, to incorporate a priori expectations about dis- 
continuity structure into the variational principle for surface reconstruction. So augmented, the 
variational principle establishes a beneficial cooperation between tlie interpolation process, which 
smoothly propagates shape information across regions, and the complementary discontinuity pro- 
cess, which delimits these regions. Thus it optimally reconstructs tlie pieccwise continuous surfaces 
and discontinuities simultaneously to achieve the best possible surface shape. 

As was mentioned in the previous section, tlie smoothness of tlie thin plate surface under 
tension is incompatible with any sudden transitions imposed by the scattered shape constraints. 
This implies that its potential energy of deformation is generally greater at what ought to be 
interpreted as surface discontinuities. Any local reduction in the continuity of die surface reduces 
the incompatibility and locally reduces potential energy. This can be seen from (8); Spt(v) 
considered as a function of (v, />, r), decreases as eitlier p[x, y) or t{x, y) are made zero over more 
of Q. This suggests that discontinuities can be discovered in die course of solving the variational 
principle, by allowing the surface to crease and fracture as needed to reduce the total energy 
below the minimum obtainable with a single smooth surfiice. The insertion of discontinuities must, 
however, incur some energy increase, otherwise /j(x, y) = everywhere would trivially minimize 
the energy. 
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Figure 10. Bending moment zero crossings (left) and detected depth discontinuities (right) for the reconstructed 
surface in Fig. 8. 



ITie variational continuity control approach to detecting discontinuities involves augmenting 
the original energy functional £pr(v) = Spt(v) + P{v) with tlie discontinuity functional P{p,t) 
(explained shortly) to obtain tlie new variational principle 



Find u, p, and r such that 



S{u,p,t) = inf e[v,p,T), 



V,p,T 



where the energy Junctional 

e{v,p,T) = $(v,p,t) -{- P{v) -\' P(p,t). 



(26) 



(27) 



The solutions u(x,y), p(x,y), and T(x,y), satisfy the diree coupled Euler-Lagrange equations, 
which express the vanishing of the first variation with respect to each independent function 

{T]Uy) ; 






6f£{u,p,f) = =f{vi, + 2vt^ + v'^^) + [1 - f]K + "y) + SpDipJ); 

Sr£{u,p,T) = =?[(kL + 2»'xj, + vly) - ("^ + "J)] + SrD{p,f). 

(28) 
Note that the first equation is identical to (9). 

The fijnctional P(p,t) maps the depth and orientation discontinuity configurations p{x,y) 
and T{x,y) into positive energies (this is analogous to the role of Spri^^) with respect to u). In its 
simplest form, the functional can increase monotonically with the totiil number of discontinuities; 
e.g., P(p, t) = J J^ /?d[l - p(x, y)] •+ /?o[l - t(x, y)] dxdy, where Pd and /?« are positive energy 
scaling parameters for the deptli and orientation discontinuity contributions, respectively. 

More interestingly, significant advantages acciiie in the detection of weak or subtle disconti- 
nuities if the functional can be designed so as to bias the solution according to generic constraints 
about the local structure of discontinuities. Usefiil constraints can, for example, be based on 
Gcstalt principles of good continuation — discontinuities tend to be arranged along contours, these 
contours tend to be continuous, etc. This may be accomplished readily by assigning potential 
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Figure 11. Rcconstmcted surfaces and disconlinuities. 



energies to various local discontinuity configurations on the (i,j) grid of nodes for tlie discrete 
problem. Encodings of local edge configurations that favor good continuation have been employed, 
for instance, in relaxation labeling curve enhancement processes [Zucker et al, 1977] and in Markov 
random field image restoration models [Geman and Geman, 1985]. 

In our current implementation, the discrete discontinuity flmctional is a weighted nodal sum 
of potential energy quanta D^j and O^^ over depdi and orientation discontinuity configurations 
respectively: 

D^p", r") = 2][/3iA^,(/) + PoOUr% (29) 

We employ a heuristic encoding which favors die formation of continuous and smoothly curving 
contours by locally assigning higher energies to isolated discontinuities, terminations, sharp bends, 
junctions, and regions. Fig. 12 illusu^atcs some of die configurations, and the (numeric) energies 
associated with them and their rotationally symmetric counterparts. Just as for computational 
molecules, the circles denote nodes («,i), while die X's denote discontinuities (positions where p^ 
or T^ are 0). The quanta in Fig. 12(a) constitute D^- for dcpdi discontinuities, which occur on 
links between nodes as explained previously. They are equivalent to die configurations found in 
[Geman and Geman, 1985]. Fig. 12(b) depicts some of Uic orientation discontinuity configurations 
encoded by Oj* •. Orientation discontinuities coincide with nodes. 
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Figure t2. Some local configurations and associated energy quanta for depth discontinuities (a) and orientation 
discontinuities (b). 



Tlie discrete variational principle simultaneously governs the values of the displacement nodal 
variables of the surface as well as the nodal variables in die discontinuity map. Although die energy 
functional Spj(v) has a unique minimum (given tJic conditions of Sec. 2.7) for fixed p and r, diis 
is no longer die case for (f (<;,/), r) which allows variation of Uic continuity control functions in the 
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minimization. The nonconvcxity of the energy landscape makes tliis a much more difficult problem 
to solve numerically. In the discontinuity detection experiments to be presented in Sec. 6.5, we 
propose a strategy for efficiently obtaining good, though not necessarily optimal solutions. 



5. Overview of the Multiresolution Surface Reconstruction Algorithm 



Application of the finite element yields die discrete problem of solving a linear system of simul- 
taneous equations. This system has computationally desirable properties; i.e., its matrix is sparse, 
banded, symmetric, as well as positive definite (for fixed p(x,y) and tau(x,y)) when the available 
constraints satisfy the conditions for a well-posed problem. The sparseness of the matrix, a direct 
consequence of the local support of die finite element, is evident from the nodal equation of an 
interior node; rows associated with interior nodes have only 13 nonzero entries, while nodes at and 
near discontinuities have even fewer. The N x N matrix however, tends to be extremely large in 
practice, since the number of pixels iV in a typical image can range from 10^ to 10^ or greater. 
This combination of properties suggests the application of iterative techniques such as (parallel) 
Jacobi or (sequential) Gauss-Seidel relaxation mediods [Hagcman and Young, 1981]. Relaxation 
methods lead to distributed algoritlims, and the parallel variants may be implemented concurrently 
on networks of many simple, locally-interconnected processors. 

5.1. Nodal Relaxation Computations 

A local-support nodal relaxation computation can be obtained at node {i,j) by expressing uf • 
in terms of the remaining variables in the nodal equation determined by the local structure of 
constraints and discontinuities. The nodal relaxation computation may be constiiictcd automatically 
by applying our simple rules governing the summation of basic computational molecules: 

(i) Plate, deptli constraint, and orientation constraint molecules sum at interior (non-discontinuity) 
nodes. 

(ii) Membrane and depth constraint molecules sum at orientation discontinuity nodes. 

(iii) Orientation discontinuities inhibit plate and orientation constraint molecules. 

(iv) Depth discontinuities inhibit all basic molecules. 

For instance, at a depth constraint node away from discontinuities, tlie Gauss-Scidel relaxation 
computation becomes 

"8 
u 



/j2 V". -2,3 + "t+2,y + "t,j- 2 + "t,i+2J + «d:,yat,j , 

where we have suppressed the discretization superscript h and instead introduced the bracketed 
iteration indices. At an unconstrained dcpdi discontinuity, we obtain 

(nf 1) _ 1 / (n + l) (n) (n+1) (n) \ , . 

"«■ J ~ 4 Vi ho + "»■ f i,j '^ "f,i-i + "t,yfi j • (31) 

Note tliat tlie nodal relaxation computations do not change from one iteration to the next, so long 
as the influencing constraints or discontinuities remain unperturbed. 
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5.2. Multiresolution Relaxation 

A serious problem with iterative techniques, in general, is their slow convergence rates for large 
problems. This inherent inefficiency is due to the fact that information must propagate incrementally 
across large representations from nodes to their near neighbors in accordance with the nodal 
relaxation fomiulas.^ We have developed highly efficient iterative algorithms that overcome this 
problem for surface reconstruction (Terzopoulos, 1982, 1983a] as well as for certain other visual 
problems [ferzopoulos, 1984]. l^hese algorithms achieve efficiency by exploiting multiresolution 
relaxation methods [Fedorenko, 1961; Brandt, 1977; Hackbusch and Trottenberg, 1982]. 

Briefly, the multiresolution surface reconstruction algorithm features (i) multiple representations 
of surface shape over a range of spatial resolutions, (ii) local, iterative (relaxation) processes 
that propagate smoothness constraints within each representational level, (iii) local coarse-to-fine 
(prolongation) processes that allow coarser representations to constrain finer ones, (iv) fine-to-coarse 
(restriction) processes that allow finer representations to constrain and improve the accuracy coarser 
ones, and (v) a multilevel coordination strategy that enables the hierarchy of representations and 
component processes to cooperate towards increasing the computational efficiency, usually by orders 
of magnitude. 

Fig. 13 depicts the structure of the algorithm schematically. In this particular case, only tliree 
levels are shown. Note the 2:1 resolution reduction between adjacent levels. Not only does this 
ratio simplify the component processes considerably, but it is also nearly optimal with regard to 
total computation to convergence (this is conveniently measured in machine independent work 
units, where a work unit is the amount of computation required for a relaxation iteration on the 
finest level) [Brandt, 1977]. The diagram illustrates the intralevel relaxation processes, as well as the 
fine-to-coarse restriction and coarse-to-fine prolongation processes that communicate between levels. 
The figure shows syndietically generated scattered orientation and depth constraints consistent with 
a hemispherical surface. The algorithm reconstmcts a dense representation of surface at three 
resolutions. The sparse information at any particular scale can be thought of as a set of constraints 
which defines a discrete surface approximation problem at that level. It is natural then to view the 
multiresolution surface reconstruction algorithm as iteratively solving a coupled hierarchy of discrete 
surface reconstruction problems.^ For a detailed description of the algoridim see frerzopoulos, 
1982, 1983a]. 



6. Experimental Analysis of the Algorithm 

ITie multiresolution visible-surface reconstruction algorithm was tested on a variety of data sets 
including synthetic data, stnictured light (laser) range data, automated stereopsis and photometric 
stereo data from natural images, and digital terrain model data. Some results are presented 
in this section (for ftjrlher details and examples, see frerzopoulos, 1984]). In all the examples 

^ It is possible to accelerate the basic relaxation methods so that fewer iterations are required. However, prac- 
tical accelerated methods such as the conjugate gradient method, successive overrelaxation, and Chebyshev 
semi-iteration use global procedures to determine the acceleration parameters. In parallel implementation, 
the greater complexity of the globally accelerated methods and, even more importantly, the communications 
costs of performing the global operations nullifies any potendal gains. ' 

^ A recursive multilevel coordination strategy was employed in the experiments described next. The recursive 
strategy activates only a single level at any one time. We have recently developed a concurrent strategy 
based on a multilevel variational principle [Terzopoulos, 1985b]. Concurrent coordination mainUiins all 
levels active simulUineously, thus achieving full parallelism. 
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Figure 13. The structure of the multiresolution surface reconstruction algorithm. Iterative relaxation processes 
operate at each level. Fine-to-coarsc ajid coarse-to-fine processes transfer information between levels. Synthetic 
orientiition and depth constraints input to Uie algorithm are shown at the top. The dense multiscale surface 
representation is output at the bottom. 
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presented, die intralcvel process was Gauss-Seidcl relaxation and Uie algorithm was started from 
zero initial approximations on all the levels. The border nodes on each grid were preset as depth 
discontinuity nodes to introduce natural boundary conditions which free the reconstructed surface 
on the boundary of n. 

6.1. Synthetic Data 

The first two examples involve randomly placed depth constraints. The left half of Fig. 14 
shows 15%-density constraints at three resolutions. These constraints were obtained by sampling 
a hemisphere whose z values were multiplied by a radial sinusoid, llie nodes outside the circular 
region occupied by the constraints were specified as depth discontinuities. The reconstaicted surface 
representation is shown on the right half of Fig. 14. In Fig. 15, the 15%-density depth constraints 
shown on the left arc samples of a stacked set of planar surfaces at tliree resolutions. In tliis 
example, depth discontinuities were placed along die circular arcs bounding tlie planes, and along 
tlie outer edges of the grids, llie reconstructed surface representation is shown on the right half of 
Fig. 15. This example indicates that discontinuities can be placed along arbitrary contours within Q 
to prevent surface shape from being degraded by unwanted smoothing over sharp depth changes. 

The next examples involve reconstructions from orientation constraints. The left half of Fig. 16 
shows in perspective a set of orientation constraints over a square region. On each of three scales, 
the region is divided into four quadrants each containing constant orientation constraints, and the 
nodes along their boundaries are preset as orientation discontinuities. The surfaces reconstructed 
by the three-level algoritiim are shown on die right half of the figure. Since absolute dcpdi cannot 
be determined solely from orientation constraints, a relative depdi reconstruction results, with die 
center of gravity of die resulting pyramidal surface resting near the x-y plane. 

The left half of Fig. 17 shows 30%-density scattered orientation constraints consistent with a 
hemispherical surface at diree resolutions. The reconstmcted surface representation is shown on 
the right. All nodes outside die hemispherical surface patch were specified as depth discontinuities. 
Again, the center of gravity of die surface rests near die x-y plane. 

The next examples demonstrate die integration of both depth and orientation constraints. The 
left half of Fig. 18 shows 15%-density depdi constraints consistent with a hemispherical surface at 
three resolutions. On die right are 15%-density orientation constraints consistent with die same 
surface. Nodes outside the surface have been specified as depth discontinuities. The reconstructed 
surface is shown in Fig. 19. Whereas in die previous example (Fig. 17) only relative depdi can be 
determined for lack of any depth constraints, in the present example the additional depth constraints 
enable the absolute depth of the surface to be determined at all points, hence the surface is "raised" 
to die correct height above die base plane. In addition, note diat (10%) uniformly distributed 
noise has been added to the constraint values. Widi die given constraint parameters, die surface is 
slightiy bumpy on die finest level. This can be reduced by decreasing die constraint parameters, in 
effect, loosening die springs of die physical model. 

6.2. Structured Light Data 

The multircsolution algorithm was applied to die reconstruction of several objects from raw 
range data supplied by a laser rangefindcr constaicted by P. Brou at MIT. The scan resolution 
in die y direction is half diat in die x direction. A four-level surface reconstruction algoridim 
was employed in die examples. The data was introduced as depth constraints at die finest level 
and transferred to die coarser levels by successive 2x2 averaging between levels. To expediendy 
segment the objects from the background, values smaller dian a diroshold were treated as depth 
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Figure 14. Reconstruction of a surface from depth constraints. (Grid dimensions: N^^ - Ny^ = 17, 
j^h, ^ ^fe, ^ 33 ]^h, ^ ^rh, ^ g5 Qj.j^ spacings: hi = 0.4. hi - 0.2, /13 = 0.1. Constraint 
parameters: ad^^ = 2.0 /hj. Computation: 24.25 worlc units.) 
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Figure 15. Reconstruction of planes with circular depth discontinuities from depth constraints. (Grid 
dimensions: N^' x Njj' = 22 x 17, N^' x Nj;^ = 43 x 33. N^' x iVjJ''' = 85 x C5. Grid spacings: 

hi = 0.4, /t2 = 0.2, /i3 = 0.1. Constraint parameters: ad^^ - 2.0/hj. Compuliition: 20.375 work units.) 
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Figure 16. Reconstruction of a pyramidal surface with orientation discontinuities from orientiUion constraints. 
(Grid dimensions: N^'^ = Nj;' = 17, N'^'^ = N^^ = 33, iV^'^ - Njj' - 65. Grid spacings: hi = 0.4, 

h2 = 0.2, /i3 = 0.1. Constraint parameters: ap' = a,^ = i.O/hj. Computation: 19.5 work units.) 
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Figure 17. Reconstniction of a hemispherical surface from scatlered orientation constraints. (Grid dimensions: 
j^hr ^ j^h, ^ yrj j^h, ^ j^h, ^ 33^ ^h, ^ j^h,-^ g5 Qrid spacings: hi = 0.4. h^ = 0.2. /ig = 0.1. 

Constraint parameters: ap' = Ug' = 4.0/hj. CompuUition: 22.125 work units.) 
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Figure 18. Depth constraints (left) and orientation conslrainls (right) consistent with a hemisphere at three 
resolutions. 
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Figure 19. Reconstruction from depth and orientation constraints in Fig. 18. (Grid dimensions: iV^* 
N^^ = 17, Nl'^ = Nj;-' --= 33, N^' = iV^'^ - 65. Grid spacings: hi = 0.4, h2 = 0.2. h^ = ( 

Constraint parameters aa'^^ - 2.0/hj, a,/ = aq^ - 4A)/hj. Computation: 17,75 work units.) 
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Figure 20. ReconsLruclion of a lightbulb from range data. (Finest grid dimensions: N^' x Ny* - 257 x 281. 
Grid spacings: hi = 0.8, h^ = 0.4, /is = 0.2, and /14 = 0.1. Constraint parameters: aj'*J = Q.2/hj. 
Computation: 9.78 work units.) 

discontinuities. Fig. 20 sliows tlie reconstructed surface of a lightbulb. 'Fhe algoritlim smoothes die 
noise in the data and reconstructs tlie missing points. 

6.3. Natural Image Data 

In this section, we apply the mullircsolution surface reconstruction algorithm to depth data orig- 
inating from natural images. The examples involve photometric stereo, and two binocular stereo 
algorithms applied to terrain stereopairs. 

6.3.1. Photometric Stereo Data 



Photometric stereo is a technique that uses multiple (usually 3) images of a scene from the same 
viewpoint, but with dilTcring illumination [Woodham, 1981]. Assuming that tlic surface material 
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Figure 21. Image of a matte while torus. 



is Icnown and tliat the viewer and lightsourccs arc far from the object, Uic mctliod determines 
the surface orientation from die image irradiance. Our surface reconstruction algoridim provides 
a noise resisuint technique for computing dcpdi from the surface orientation data provided by 
photometric stereo. We demonstrate this with tlie image of torus in Fig. 21. The photometric stereo 
data (generated by a system implemented at MIT by K. Ikeuchi) was introduced as orientation 
constraints on a two-level algoritlim. Aside from sporadic missing data, the constraints on die 
coarse level arc dense, whereas only every other node on die fine level is a constraint. Fig. 22 
shows the orientation data and die reconstructed torus. 

Our method for reconstructing surfaces from scattered orientation constraints can be compared 
to a variational scheme for obtaining relative depth from dense surface gradient infonnation reported 
by Horn and Brooks [1985]. Their proposed least squares integral / /(vi - p)^ + (vy - q)^dxdy 
will be recognized as being a continuous version of die orientation constraint penalty functional. By 
virtue of the additional smoothness functional Sprlv), however, our surface reconstmction algoriUim 
can deal with orientaUon constraints that are scattered. It also can integrate depdi constraints from 
other sources to arrive at absolute surface depdi. 

6.3.2. Correlation Based Stereo Data 

At die top of Fig. 23 is a stcrcopair on which Kass's [1983] correlation based stereo algoridim 
was run. The output of the stereo algoridim is shown on the lower left, with brightness proportional 
to disparity. 'I he algoridim has failed to produce a match in the neutral grey patches, so disparity 
is' unknown in dicse areas. To apply the mukircsolution algoridim, the disparity data on the finest 
level were reduced by factors of two, dirough averaging, to dirce coarser levels. RelaUvely small 
constraint parameter values were clioscn in order to counteract die potentially detrimental effects 
of false matches and noise in die disparity data. The rcconstaictions on die diree coarsest levels 
are shown as 3D plots in Fig. 24 (the finest level was too dense to represent diis way). Fig. 25 
shows isoelevauon contour maps of the solution on all levels. 
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Figure 22. Reconstruction of the torus (right) from the orienUition constraints provided by photometric stereo 
(left). (Grid dimensions iV^ = iV^*> = 51 and N^^ = N^^ = 101. Constraint parameters: a''^ = A.O/hj. 
Computation: 52.0 work units.) 

6.3.3. Feature Based Stereo Data 

The next example involves disparity constraints generated by the MPG stereo algorithm 
[Crimson, 1985]. A three-channel version of the stereo algoritlim was run on the stcreopair at 
the top of Fig. 26. The output of tlie stereo algorithm is shown on tlic lower part of the figure. 
Disparity information is provided only along zero crossing contours at the three finest scales. In 
tlie figure, the darkness along contours is proportional to disparity. This disparity data was input to 
a four-level surface reconstruction algorithm. The constraints on the coarsest level were derived by 
averaging tlie constraints from the next finer level. 'The reconstructions on tlie three coarsest levels 
are shown as 3D plots in b'ig. 27. Fig. 28 shows isoelevation contour maps of the solution on all 
levels. 



6.4. Digital Terrain Map Data 
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Figure 23. Nalural terrain stereopair (top) and output of Kass' stereo algorithm (bottom). The images were 
25G X 25C pixels, quantized to 256 levels (provided by the US Defense Mapping Agency). 



Terzopoulos 



Computing Visible-Surface Representations 



39 




Figure 24. Reconstruction of terrain in Fig. 23. (Grid dimensions: N^' =^ N^' ^ 33. iV^^ = N^-* = 65, 
^hs ^ j^h, ^ 129, iV^ = N^' = 257. Grid spacings: hi = 0.8, /i2 = 0.4, hz = 0.2, h^ = 0.1. 
Constraint parameters: a^'*-' - 0.0J//iy. Compulation: 29.0 work units.) 
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Figure 25. Isoelevation contour maps of Uie reconstructed terrain in Fig. 24. 



A four-level surface, reconstruction algorithm was applied to contoured terrain elevation data. 
A contour map of tlie iilack River Gorges (published by tlie UK Ministry of Defense) was 
digitized manually on a digitizing tablet by J. Mahoney. The 256 x 256 digital contour array is 
shown at the top of Fig. 29. The constraints input to the algoritlim are shown at the bottom 
of the figure. The elevation of the contours is proportional to brightness. Local averaging was 
used to derive the constraints on the coarser grids from those on the finest grid. The terrain 
reconstructions on die three coarsest levels are shown as 3D plots in Fig. 30 (the finest level is too 
dense to represent this way). Fig. 31 shows isoelevation contour plots of the reconstructed terrain 
on all levels. The reconstructed contours on the finest level can be compared subjectively with 
tlie digitized contours in Fig. 29, but note the reconstructed contours depict elevations half way 
between tlie original constraint contours for an unbiased comparison. The reconstructed contours 
are somewhat smoother than tlie (prcdigitized) contours in the original map — the jaggcdness 
introduced by manual digitization has been reduced. The extent of the smoothing can be regulated 
by adjusting the constraint parameters. Shaded image renditions of die reconstructed terrain using 
reflectance map techniques for hill shading [Horn, 1981] are shown at the bottom of F'ig. 31. Terrain 
reconstructions using the thin plate surface under tension model were compared to reconstructions 
using the simpler membrane spline model (Faplacian smoothing). The former gives good results, 
whereas the latter generally suffers from insufficient smoothness and produces flat spots across 
terrain peaks [Ferzopoulos, 1984] (see also [iJolondi el al, 1976]). 
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Figure 26. Natural terrain stereopair (top) and output of the MPG stereo algorithm (bottom). The images 
were 512 x 512 pixels, quantized to 256 levels (provided by the US Army Engineer Topographic Labs). 
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Figure 27. Reconstruction of data in Fig. 26 on the three coarsest levels. (Grid dimensions: N'^^ =- N^^ = 33, 
^fe, = Nj^^ = 65. N'^'^ = N'y'^ = 129, iV^ = N^' = 257. Grid spacings: hi = 0.8, /iz =^ 0.4. /13 = 0.2, 
/14 = 0.1. Constraint parameters: aj^^ - Q.Ol/h'j. Computation: 31.0 work units.) 



Terzopoulos 



Computing Visible-Surface Representations 



43 








.•^^-•r^i 



LlW 



i^^W, 





Figure 28. Isoelevation contour maps of the reconstructed terrain in Fig. 27. 

6.5. Discontinuity Detection Experiments 

The foregoing examples have shown that the surface reconstniction algorithm can handle disconti- 
nuities that are prcspecified. In tlie next two sections, wc present examples involving tlie automatic 
detection of discontinuities. 

6.5.1. The Rcgularization Approach 



The aerial view stereopair of Fig. 32 was input to tlie MPG stereo algorithm which generated 
the sparse disparity map shown on the top of Fig. 33. The finest level dense disparity map generated 
by a four-level surface reconstruction algorithm is shown at the lower left of the figure. Darkness 
is proportional to disparity. ITie discontinuities found from this disparity map, using a disparity 
limit Gij > id = 1 arc shown at the lower right as white contours After the detected points are 
added to the discontinuity map, the surface reconstruction algoritlim continues iterating from the 
tentative approximation on tlie left. The amount of additional computation required is relatively 
small, since tlie tentative surface is a fairly good approximation in most places. At convergence, the 
reconstructed surface has fractured along the contours to give tlie solution on the right. Portions 
of tlie main discontinuities around tlie buildings have been found, but contours are broken and 
shifted. 

The next example involves the synthesized random dot stereogram in Fig. 7. The depth 
constraints generated by a three-channel version of tlie MPG stereo algoritlim are shown in l-ig. 34 
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Figure 29. Digitized contour data (top) and constrain ts (bottom). Tlie patcii to the lower riglit represents a 
lake. 
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Figure 30. Reconslruclion of data in Fig. 29. The terrain reconstruction on the lliree coarsest levels is 
represented as 3D surface plots. (Grid Dimensions: A^^ -^ A^^'" =^ 33, Nl'' = Nj^'' = C5, N^' = N^^ ---- 
129, N^* = N^* = 257. Grid spacings: hi = 0.8, h-i ^: 0.4. /13 = 0.2, /14 = 0.1. Constraint parameters: 
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Figure 31. Isoelcvation contour map (lop) of the data in Fig. 30 and shaded representations of the reconstructed 
terrain (bottom). 
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Figure 32. Aerial view of a hospital complex. The stereopair was provided by the UBC Faculty of Forestry. 
Images are 320 x 320 pixels. 

(the finest level dimensions are 320 x 320). The constraints on tlie coarsest level were obtained by 
averaging those on the next finer level. Fig. 35 shows the smooth disparity maps initially computed 
by a four-level surface reconstruction algorithm. Fig. 36 shows the discontinuities detected from 
these maps witli tj = 1. The discontinuities have been superimposed onto the final disparity maps 
in Fig. 37. Better performance is observed in tliis case due to the simpler surface structure, but the 
contours, while mostly intact, are quite ragged. 

In general, not detecting true discontinuities affects surface shape more adversely over larger 
regions tlian introducing some spurious ones within a continuous surfiice. Discontinuity points 
are missed by llie thresholding operation, and no adjustment of the global limit can be expected 
to produce perfect results. Note, however, that the surface reconstruction algoritlim does not 
break down. Ratlier, the reconstructed surface degrades as it "leaks" through the gaps. The 
discontinuity detection procedure ma(y be improved by allowing the disparity limit to vary spatially, 
or by modifying it during multiple passes. On the first pass, surface shape is poorest, so a fairiy 
conservative limit should be set to reduce tlie number of false detections. Conservative limits fail 
to detect many discontinuities, but as more discontinuities are identified, surface shape improves 
and limits can be lowered in subsequent passes to find the less prominent discontinuities. 

6.5.2. The Variational Continuity Control Approach 



A multipass scheme is also employed in the variational continuity control approach to efficiently 
obtain good solutions. An example will be used to explain the strategy. Fig. 38 shows depth 
constraints randomly sampled from a set of sloping planes that forni discontinuities along their 
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Figure 33. Discontinuities in the aerial stereogram. Disparity contours generated by stereo algorithm (top), 
full disparity map generated by the surface reconstruction algorithm at the finest level (lower left), and detected 
discontinuities superimposed on the disparity map (lower right). 
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Figure 34. Depth constraints for the random dot stereogram. 



extremities. A single continuous surface can be reconstructed from tliese constraints as shown, but 
it smooths over tlic depth discontinuities and rounds out Uic orientation discontinuities. 

The algorithm finds both kinds of discontinuities and reconstaicts a surfiice which preserves 
them. When the surface smooths through a deptli discontinuity, two spurious regions of high 
curvature border tlie discontinuity. These spurious regions can easily be mistaken for orientation 
discontinuities, i b avoid tliis unwanted interaction which can substantially slow down the optimiza- 
tion process, tlie algoritlim postpones the orientation discontinuity detection phase until all depth 
discontinuities have been found. The surface evolves in several steps over which the parameters 
fi^^ and p^ in (29) are modified. Each step consists of first flipping the value of the continuity 
control parameter (/?J* • or r/' ) from to 1 or conversely, if this lowers the energy (27), and then 
nmning the rcconstmction algorithm to convergence (which always results in equilibi"ium, since tlie 
variational principle is convex for fixed p^- and t^X 

For depth discontinuities, /?J is initially set to a high value that heavily penalizes their insertion, 
tlien lowered in steps. This strategy of least commitment finds the prominent discontinuities earliest, 
improving the sinface as it does so, and leaves tlie more subtle ones for later. It results in the 
flipping of relatively few variables in each stage, hence the solution is obtained efliciently. Beginning 
with the continuous surface V\%. 39(a), Fig. 39(b-d) illustrates the steps of tlie evolving discontinuity 
detection process, during which discontinuities are determined with increasing accuracy as ft^ is 
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Figure 35. Full disparity maps without discontinuities. 



lowered. The energy can be lowered further still if p^ is then increased slightly to eliminate 
spurious discontinuities in Fig. 39(d). Note that since the surfaces have now separated, a very large 
increase would be needed to flip a true discontinuity point (a hysteresis effect). The improved 
surface in Fig. 39(e) results. Next, tlie orientation discontinuity detection phase is activated and it 
nms in tlie same way, but modifies ft^. In this example, tlie orientation discontinuities are found 
in only one step. 

Fig. 39(f) shows tlie final solution. The depth and orientation discontinuities have been made 
explicit and are preserved by the reconstructed surface. Incidentally, tlie global optimum of the 
variational principle has been found in tliis example; however, this procedure can generally be 
expected to yield good, though not necessarily optimal approximations. Its main attractions are 
that it is deterministic and efficient. 



7. Discussion and Research Directions 



Several issues concerning tlie framework for computing visible-surface representations are discussed 
in tliis section, and directions for future research are suggested. The discussion focuses on 
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Figure 36. Detected discontinuities. 



discontinuity detection, choosing constraint parameters, iiandling rivalries in constraints, grouping 
constraints, invariancc properties of the surface reconstaiction model, and visible-surface analysis, 
frerzopoulos, 1984, Ch. 11] contains a more extensive treatment of tliese and other issues, including 
multircsolution relative depth representations of surfaces, and the possibility of computing visible- 
surface representations "instantaneously" by analog networks. 



7.1. On Discontinuity Detection 



Some recent work in image restoration is of relevance to the problem of piecewise continuous 
surface reconstruction. A piecewise. constant image model employed by Blake [1983] for image 
reconstruction is interesting in tliat it incoiporatcs "weak constraints" which can be broken at a cost. 
The resulting optimization problem is related to our variational continuity control approach, but 
more restricted. Blake used an adaptive method, which he referred to as "graduated nonconvexity," 
to obtain good solutions to the nonconvex problem. It has not been established however whether 
this interesting method applies to tlic sparse data case as well. 

Geman and Geman [1985] used Markov random field models with associated Gibbs distributions 
to restore piecewise constant images corrupted by additive Gaussian noise. The restoration seeks 
a maximum a posteriori estimate of the original image, given the degraded image, and includes 
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Figure 37. Disconlinuities aiid final disparity maps. 



an explicit "line process" that estimates the locations of step edges in intensity. This work was 
restricted to dense image data. The Gcmans' approach was adopted with encouraging results to 
surface reconstruction from sparse depth constraints by Marroquin [1984]. His markov random field 
model, while less restrictive than tlic Gcmans' pieccwise constant one, in fact models a membrane 
spline whose smoothness is insufficient for computing visible-surfoce representations. A line process 
essentially equivalent to the Gcmans' was incorporated to estimate deptli discontinuities, llie 
numerical solution strategy in both of the above studies was stochastic optimization using the 
Metropolis algorithm and simulated annealing to optimize the nonconvex functional [Kirkpatrick 
el al., 1983]. This strategy can find optimal solutions, but for such large reconstruction problems it 
has been observed to converge notoriously slowly. Based on our experience, we believe that it can 
be accelerated, perhaps enough to make it practical, through ilie use of multiresolution processing. 

Obviously, the line processes used in the above work as well as our own encoding of disconti- 
nuity contour configurations is unpleasingly heuristic and in need of refinement. The discontinuity 
map can be augmented by nodal variables to encode the local orientations of the curvilinear ele- 
ments to a higher degree of accuracy. Such an encoding is employed by Zuckcr and Parent [1984] 
in an optimization (relaxation labeling) approach to finding contours in images. It appears that 
ideas from their work can also be applied to finding surface discontinuities within our framework. 

A promising possibility is to employ ID controllcd-continuity stabilizers as formal models 
of smoothness constraints along surface discontinuity contours in the x-y plane. A ftinctional 
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Figure 38. Scattered depth constraints consistent with sloping planes meeting discontinuously (top) and the 
smooth reconstructed surface (bottom). 

that naturally comes to mind is the curvilinear analog of die thin plate surface under tension: 
f^ /3i, (s) {fta(s) {d^c/ds^) + [l-fia (•'<)1 (dc/ds) } (is, where s denotes arc length along discontinuity 
contours c E C. Here, pi allows breaks, while pa allows angles (tanjent discontinuities) to form 
in the discontinuity contours. Again, additional energy penalties must be associated with these 
occurrences. Given our finite element representation of surfaces, curvilinear finite elements are 
tlie natural local representation for discontinuity contours. The combined variational principle has 
both a surface component and an analogous contour component. Although technically nontrivial, 
a fonnulation of surface rcconstniction generalized along these lines has very strong appeal. 

7.2. On Constraints —• Parameters, Rivalries, and Grouping 



The constraint (spring) parameters ofl'er tlic flexibility to individually tune the coerciveness of each 
constraint on the reconstructed surface. In the special case of Gaussian error distributions, the 
parameters should be inversely proportional to the expected variances (a,- = l/Xoj). It ought to be 
possible for the low-level visual processes to associate a variance estimate or confidence with each 
constraint that they provide. In general, however, it's not obvious how to choose tlie constraint 
parameters optimally. 

The constant of proportionality A~^ can also be used to tune the overall smoothness of 
die reconstructed surface. Cross validation techniques may be used to set A optimally (e.g., 
[Wahba and Wendelbcrger, 1980]). The basic criterion is to choose A so as to minimize over all 
constraints tlic (weighted) discrepancy between each constraint and its value as estimated from the 
surface reconstructed using tlie remaining constraints. Unfortunately, this involves computationally 
expensive sequential algorithms. Interestingly, the continuous tuning of surface smoothness is 
analogous to the scale space filtering technique proposed by Witkin [1983] with the added attraction 
that it can be applied to scattered data. 

Although the variational principle was designed to account for measurement errors in the 
constraints, the possibility of massive rivalries between constraints from different sources, such as 
stereopsis and analysis of motion processes, was disregarded. Massive rivalries are unnatural visual 
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Figure 39. Evolution of the discontinuity detection process. 



phenomena that can nevertheless occur, especially under contrived conditions, and they often lead 
to multistable percepts [Attneave, 1971]. The framework can potentially accommodate rivalries with 
a mechanism that inhibits individual or entire sets of constraints by nullifying selected constraint 
parameters. This mechanism can be activated by a global arbitrator which monitors the contents 
of visible-surface representations to detect rivalries may also have access to higher level knowledge 
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about the scene. ITie arbitrator's influence can account for multistability. 

A particular type of rivalry arises from transparent surfaces. For instance, a surface such as a 
dirty window in front of a background scene would lead to two well defined populations of depth 
(and orientation) constraints over the same visual angle, one from the window, the other from the 
background. A transparency interpretation can be arrived at by an arbitrator which monitors the 
surface reconstruction process looking for high approximation error between surface and constraints 
over a significant area. Under these conditions the arbitrator can trigger a constraint grouping 
process which clusters the constraints into two populations, based on depth values, say. Multiple 
surfaces can then be reconstructed over the same visual area for each resulting constraint population. 
This scheme has been applied on a transparent surface random dot stereogram (Terzopoulos, 1984, 
Ch. 11]. 

7.3. On Invariancc Properties of the Surface Model 



As a transformation from sparse constraints to dense surfaces, the thin plate under tension model 
can be shown to be invariant under (i.e. tommutes with) certain image plane transformations 
applied to the constraints; namely, translations, rotations, and similarity transformations. This 
implies that surface shapes will be preserved through rigid motions of the scene or viewpoint 
parallel to the image plane or along the view direction. These are essential invariance properties 
for visible-surface reconstruction [Terzopoulos, 1982]. 

Note, however, that the thin plate spline, characterized by the small deflection approximation 
/ / ^iz + 2viy + Vyj, dx dy to tlie bending energy density of a thin plate, is not invariant under 
arbitrary 3D tiansformations of the constraints. Thus, surface interpolation using this expression is 
not invariant under changes in the view direction, as Blake [1984] points out. He shows that rotating 
tlic view direction induces the ID analog of the thin plate spline to "wobble," and he demonstrates 
tliat this effect is most pronounced as the (continuous) spline is inclined sharply with respect to the 
viewer or is forced to bend sharply. Blake views this as a problem tliat should be eliminated by 
employing the large deflection bending energy of the thin plate, a convex combination of the mean 
and Gaussian curvatures of the surface v(x,y), which is view direction invariant. 

Although Spr(v) can also be made view direction invariant by employing the large deflection 
counterparts for tlie tliin plate and membrane bending energies, this approach has a serious 
technical drawback [Terzopoulos, 1984]: The large deflection formulas lead to an extremely difficult 
nonlinear problem (e.g., the large deflection equations for the thin plate are two coupled nonlinear 
fourtli-order partial differential equations known as Von Karmann's equations [Szilard, 1974]). 

Fortunately, the surface rcconstniction model, as it stands, is not hampered by die lack of 
view direction invariance because the available constraints are usually sufficiently dense in practice 
to tightly determine surface shape; as the view direction is varied, tlie reconstructed surface would 
vary negligibly (note that Blake's experiments reveal a significant wobble eff^ect just in tlie case 
of extremely sparse constraints). Inirtliennoro, the explicit introduction of depth and orientation 
discontinuities alleviates much of the wobble precisely at those places where Blake's experiments 
show it to be most pronounced on a globally continuous surface. An interesting psychophysical 
experiment would be to determine whether there might be some slight variance in the surfaces 
perceived by humans viewing sparse random dot stereograms while the dots undergo simulated 
rigid 3D transformations and, if so, whetlier tlie variations arc consistent with tlie reconstruction 
model (J. Mayhew, personal communication). 
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7.4. On Visible-Surface Analysis 

The visible-surface representation is an intemiediate and volatile description of the 3D surfaces 
in scenes. It drives ensuing processes which generate stable higher-level representations of shape 
that are better tuned to object recognition. The processing begins with visible-surface analysis 
whose goal is to abstract from the numeric, viewer-centered representation a rich set of more 
symbolic, object-centered features that are stable through viewpoint changes. The extraction of 
geometric surface features is facilitated by the dense shape information provided by visible-surface 
representations. 

A promising approach to visible-surface analysis is to apply concepts from differential geometry 
[do Canno, 1976]. For instance, a surface's intrinsic geometry (including Gaussian curvature, 
geodesies, etc.) is determined completely by the first fundamental form, which defines arc length 
over die surface. It's extrinsic geometry (including normal curvature, principal curvatures, etc.) are 
determined by tlie second Hmdamental form, which describes the deviation of the surface from the 
local tangent plane. ITic fundamental dicorcm of the local dicory of surfaces (usually attributed 
to Bonnet) states that Uie analytic study of surface properties consists of the study of tlie two 
fundamental forms; i.e., the six fundamental tensor coefficients (which arc not all independent) as 
functions of die two independent parameters of die surface. The fimdamental forms arc invariant 
under changes in the parameterization, and together they determine surface shape up to rigid body 
transformations. These properties make Uiem ideal foundations for object centered symbolic surface 
representadons. 

The visible-surface representation makes it possible to estimate the first and second funda- 
mental forms on a point-by-point basis over die entire visible surface. .The finite element shape 
representation reduces die computation of crucial local surface features such as die Gaussian curva- 
ture, principal curvatures, and principal directions to the evaluation of simple algebraic expressions 
of neighboring nodal variables (see [Terzopoulos, 1984, Ch. 11] for derivations). It is then a simple 
step to detennine die elliptic, hyperbolic, parabolic, umbilic, and planar points, as well as geodesies, 
asymptotes, and lines of curvature. 

For example, Fig. 39 shows the reconstructed surface of a lightbulb. Fig. 40 shows the Gaussian 
curvature K{x,y) computed for die reconstaicted lightbulb surface of Fig. 20. llie elliptic points 
{K > 0) are shown in white, die hyperbolic points (K < 0) are shown in black, and the parabolic 
{K — 0) points separate die two regions. Note the alternation in the sign of curvature at die 
screw mount. Fig. 41 plots die computed field of principal directions for the lightbulb at the 
two coarsest scales. These demonstrations illustrate die feasibility of reliably computing from diese 
representations higher-order intrinsic and extrinsic properties of surface shape. The reliability can 
be attributed to the regularizing properties of the thin plate surfiice under tension which overcomes 
die potentially detrimental effx^cts of noise in the data, while preserving discontinuities. For further 
analysis of die kinds of features diat can be computed from dense, numeric, representations of 
surfaces see, e.g., [Brady el a!., 1985] or [Medioni and Nevada, 1984]. 



8. Conclusion 



Constraints on surface shape, contributed by multiple low-level visual processes, can be computed 
reliably at multiple resolutions, but only at scattered locations in the field of view. Subsequent 
visual processing can be facilitated substantially if the scattered constraints arc transformed into 
visible-surface representations diat make surface shape explicit everywhere. 'I'o accomplish this 
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Figure 40. Elliptic (white) and hyperbolic (black) points of the reconstructed lightbulb at four scales. 



effectively, information must be integrated over multiple visual modalities and fijsed across multiple 
scales of resolution. 

In this paper, we have developed a computational theory of visible-surface representations. 
Within a unified computational framework, formal solutions were offered to ftindamental problems 
of reconstructing visible surfaces: (i) integrating constraints on die depth and orientation of surfaces 
across various modalities and scales, (ii) interpolating surface shape infonnation into (piccewise) 
smooth surfaces, (iii) discovering discontinuities in surface depdi and orientation and enabling them 
to restrict interpolation, and (iv) efficiently maintaining consistency in distributed, multiresolution 
visible-surface represcntadons. 

A visible-surface reconstruction algoridim implements the framework. Extensive testing has 
shown it to be viable. 1'he algorithm coordinates cooperative processes within a multiresolution 
hierarchy of surface representations to dramatically increase computational efficiency. It is well 
suited to implementation on massively parallel networks of simple, locally interconnected processors. 
Such computational networks are suggestive of biological mechanisms and are also well suited to 
VLSI technology. 
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Figure 41. Principal directions for tlie reconstructed lightbulb at the two coarsest scales. The directions of 
greatest curvature are shown on ihe led Those of least curvature are shown on the right 
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